#include "include/ODESolver.h"
#include <ctime>

using namespace std;

#define EXP 2.718281828

double lambda = -100;
double eta = 1.5;
// double k = 0.01;
// double T = 1;

// // 需要注意 N = T / k 的计算方法可能不准
// int N = 100;

// 应该这样计算
double T = 1;
int N = 1000;
double k = T / N;

Vector f(Vector x, Real t)
{
    return lambda * (x[0] - cos(t)) - sin(t);
}

Real u(Real t)
{
    return pow(EXP, lambda * t) * (eta - 1) + cos(t);
}

/**
 * @file main.cpp
 * @author xingyifan
 * @date 2022-07-27 22:07
 *
 * @description: 经典 RK 方法
 */

void testClassicalRK()
{
    cout << "Testing classical RK" << endl;
    ODESolver solver;
    ClassicalRK method(k);
    // 注意要多给一个
    TimeGrid G(N + 1, 0);
    G[0] = {eta};
    solver.solve(f, G, &method);
    cout << G.back() << endl;
    cout << u(T) << endl;
    cout << method.t_ << endl;
}

/**
 * @file main.cpp
 * @author xingyifan
 * @date 2022-07-27 22:09
 *
 * @description: ESDIRK 方法
 */

void testESDIRK()
{
    cout << "Testing ESDIRK" << endl;
    ODESolver solver;
    ESDIRK method(k);
    TimeGrid G(N + 1, 0);
    G[0] = {eta};
    solver.solve(f, G, &method);
    cout << G.back() << endl;
    cout << u(T) << endl;
}

/**
 * @file main.cpp
 * @author xingyifan
 * @date 2022-07-27 22:10
 *
 * @description: GaussLegendre 方法
 */

void testGaussLegendre()
{
    cout << "Testing GaussLegendre" << endl;
    GaussLegendre method(3, k);
    ODESolver solver;
    TimeGrid G(N + 1, 0);
    G[0] = {eta};
    solver.solve(f, G, &method);
    cout << setprecision(10) << G.back() << endl;
    cout << setprecision(10) << u(T) << endl;
}

int main(int argc, char **argv)
{
    testClassicalRK();
    testESDIRK();
    testGaussLegendre();

    return 0;
}